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Using numerical techniques, we study the collapse of a scalar field configuration in the Newtonian 
limit of the spherically symmetric Einstein-Klein-Gordon (EKG) system, which results in the so 
called Schrodinger-Newton (SN) set of equations. We present the numerical code developed to evolve 
the SN system and topics related, like equilibrium configurations and boundary conditions. Also, 
we analyze the evolution of different initial configurations and the physical quantities associated to 
them. In particular, we readdress the issue of the gravitational cooling mechanism for Newtonian 
systems and find that all systems settle down onto a 0-node equilibrium configuration. 
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I. INTRODUCTION 

In a previous paper of oursU, we studied the forma- 
tion of a gravitationally bounded object comprised of 
scalar particles, under the influence of Newtonian gravity. 
The dynamics of the system is described by the coupled 
Schrodinger-Newton (SN) system of equations, which is 
nothing but the weak field limit of its general relativistic 
counterpart, the Einstein Klein Gordon (EKG) system. 

As at the moment we have no hints to finding an an- 
alytic solution for the evolution, we found necessary to 
develop numerical techniques to study the formation pro- 
cess of the scalar objects. The study of the dynami- 
cal properties of the fully time-dependent SN system has 
been done before in the literatureHHIiQ, 

but more is 

needed in order to understand the gravitational collapse 
of a weakly gravitating scalar field. 

The main aim of this paper is to perform an exhaustive 
numerical study of the collapse and evolution of a spher- 
ically symmetric scalar object in the Newtonian regime. 
Here, we develop a numerical strategy to evolve the SN 
system, and study important issues like the stability and 
the formation process of gravitationally bound scalar sys- 
tems, a topic that has recently become attractive in Cos- 
mology Ef I Hi 

A summary of the paper is as follows. In Sec. [HI we 
present the relativistic EKG and Newtonian SN equa- 
tions that describe the evolution of a self-gravitating 
scalar field in the spherically symmetric case. Corre- 
spondingly, it is described how the EKG and the SN 



equations are solved in order to obtain regular and 
asymptotically flat solutions. These solutions arc known 
as boson stars and oscillatons, respectively, for complex 
and real scalar fields. However, we focus our attention 
on their corresponding weak field limit, SN system and 
its properties. 

In Sec. IIIII we present an appropriate numerical code 
to evolve the SN system, providing details about the algo- 
rithms used. The issues concerning the physical bound- 
ary conditions imposed on the SN system and the accu- 
racy of the code are of particular importance. 

The results of the numerical evolutions are given in 
Sec . 11 VI We systematically test the boundary conditions, 
the convergence of the numerical solutions and how the 
latter reproduces the expected results for the equilibrium 
configurations of the SN system. 

Sec.[V]is devoted to the study of the gravitational cool- 
ing mechanism, first described inQ. Finally, the conclu- 
sions are collected in Sec. I VII 
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II. MATHEMATICAL BACKGROUND 

Here we describe the mathematical basis to deal with 
classical scalar fields, complex and real, in both the rel- 
ativistic and Newtonian limits. The latter is given in 
much more detail as it will be our system of interest for 
the rest of this paper. 

To begin with, we write the equations that describe a 
scalar system within General Relativity, which are the 
coupled Einstein-Klein-Gordon equations (EKG). For 
simplicity, we deal only with the spherically symmetric 
case for a single scalar fluctuation. Then, we describe 
how the EKG equations can be solved to obtain regu- 
lar and asymptotically flat solutions, which are known in 
the literature as boson stars and oscillatons. Some com- 
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mcnts on the stability of these relativistic solutions are 
also given. 

After that, we obtain the Newtonian limit of the EKG 
system through a post-Newtonian procedure, which 
yields the so called Schrodinger-Newton (SN) equations. 
The Newtonian limit is quite interesting by itself be- 
cause many physical quantities can be defined and mea- 
sured, quantities that are useful to describe the system 
in a detailed manner. Moreover, the SN system obeys a 
scale transformation such that the study of all the pos- 
sible equilibrium configurations is reduced to the study 
of properly sized profiles, and this includes the evolution 
process too. 

The SN system can be solved to find stationary so- 
lutions that we shall call Newtonian equilibrium config- 
urations. These arc called Newtonian boson stars and 
Newtonian oscillatons, following the relativistic nomen- 
clature. Though Newtonian oscillatons are described by 
a larger set of equations than Newtonian bosons stars, 
their dynamical evolution is governed by the same SN 
system. Finally, we will also calculate the perturbations 
of the 0-node equilibrium configurations. These pertur- 
bations are a particular imprint of 0-node solutions of 
the SN system, and also indicate that the latter are in- 
trinsically stable. This information will be useful for the 
numerical studies performed in the following sections. 



A. The Einstein— Klein— Gordon equations 

The energy-momentum tensor of a complex scalar 
field $( c ) endowed with a scalar potential l/(|$( c )|) = 
m 2 |$( c )| 2 , reads 



rp( C ) 



<T>( C )<T)( C )* + cf>( c )*<J>( c 



$(<0><r$(c)* 



l$( c )l 2 



(1) 



whereas for a real scalar field $( r ' endowed with a poten- 
tial = (1/2)to 2 $M 2 , is given by 

T$ = *<*><l>« - (<$>^^) + ™ 2 $« 2 ) . (2) 

The parameter m is interpreted in both cases as the mass 
of the scalar particles. 

As we will work on the Newtonian limit of the scalar 
fields, we represent them using the same variables since 
much of the analytical treatment is similar for both type 
of fields. It is, however, easy to distinguish each case, 
as we will refer to complex fields wherever a complex 
conjugation appears in the formulas. 

The equation of motion governing the evolution of the 
scalar field is the Klein-Gordon (KG) equation, which 
appears from the conservation of the energy-momentum 
tensor 

T (cr)^ = Hcncej the KG equation is written, 
respectively for the complex and real cases, as 



where □ = -^==e^ [\/— gg^ v d v ] is the covariant 
d'Alambertian operator. 

Within General Relativity, the evolution of the sys- 
tem is governed by the KG equation coupled to the Ein- 
stein equations. For simplicity, we consider the spheri- 
cally symmetric case in the polar-areal slicing, for which 
the metric is written in the form 



ds 2 



-a 2 dt 2 + a 2 dr 2 



(4) 



□ $(c,r) _ m 2 $ (c,r) = ^ 



(3) 



where a(t, r) and a(t, r) are functions to be determined 
sclf-consistcntly from the matter distribution. The Ein- 
stein equations are written as usually, G^ v = 8nGT^, 
being G^ the Einstein tensor corresponding to the met- 
ric Eq. (0J. The term on the right hand side is the 
energy-momentum tensor Q or J2J. Notice that we are 
using units in which c = h = 1, and then the Planck mass 
is given by mpi = G -1 / 2 . 



B. Relativistic equilibrium configurations 

Non-trivial solutions of the (spherically symmetric) 
EKG system exist that are regular everywhere, asymp- 
totically flat, and for which the scalar field is confined to 
a finite region. These solutions are equilibrium configu- 
rations known as boson stars and oscillatons, whether we 
speak of complex or real scalar fields, respectively. We 
give here a brief description of both solutions. 



1. Boson stars 

Boson stars are the simplest solutions since the Klein- 
Gordon field admits an stationary solution of the form 
v 7 47rG &( c >(t, r) = e -lw< 4>{r) , where 4>{r) is a real function 
and uj is called the fundamental frequency of the boson 
star. Then, the metric functions are time-independent, 
which can be seen from the cancellation of the time de- 
pendence in the expression for the energy-momentum 
tensor 

The EKG system becomes a set of coupled ordinary 
differential equations that has to be solved numerically. 
By imposing the conditions of regularity at the center, 
a(r = 0) = 1 and <fi' (r = 0) = 0, together with the con- 
dition of asymptotic flatness, the EKG system becomes 
an eigenvalue problem. For each value of the field at the 
center, 4>{r = 0) = </>o, there are unique (eigen)values of 
uj and a(r = 0) = ao for which the conditions mentioned 
above are satisfied. 

The resulting boson stars have been widel y st udied in 
the literature ElllISIIlEHQHUlEE 
Hi! l20l | , including the cases of non-spherically symmetric 
equilibrium configurations 0, 12^, HjJ 0] . Their simple 
properties have also inspired some models of dark matter 
in galaxies. We will not comment more on boson stars, 
but the interested reader will find useful information in 
the references given above and in the reviews [2f|, [2(| . 
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2. Oscillatons 

For the case of a real scalar field, it is not possible to 
propose regular and stationary solutions such that the 
metric functions arc time independent. Rather, it has to 
be taken into account that all fields are time-dependent. 
For this, it is usually assumed that a Fourier expansion 
of the field and the metric coefficients suffices to se t up 
the behavior of the equilibrium configurations^. I2H l28j . 
as described next. 

First of all, we take the following Fourier expansions^ 

m 

00 

V8ttG* w (*.»•) = 5^-(r)cos(jwt), (5a) 

OO 

A(t,r) = J2 A j(r)cos(ju>t), (5b) 

OO 

C(t,r) = ]TQ(r) cos(j^), (5c) 
3=0 

and substitute them into the EKG equations. Here, 
A(t,r) = a 2 (t,r), C(t,r) = a 2 (t,r)/a 2 (t,r), and uj is 
the fundamental frequency of the system. 

A set of coupled ordinary differential equations is ob- 
tained by setting each collected Fourier coefficient to 
zero. The non-linearity of the Einstein equations shows 
up because the different modes are mixed, and then the 
equations for all the modes have to be solved simulta- 
neously. In order to avoid an infinite set of equations a 
cut-off mode is introduced, and thus expansions J5J are 
taken up to a maximum value j = j max , and then modes 
for which j > j max are eliminated by hand. 

As for the boson star case, we impose the condition of 
regularity at the center, which now reads <t>'j(r = 0) = 
and ao(r = 0) = 1 ,a,j>x(r = 0) = 0. Again, we assume 
asymptotic flatness which converts the solution of the 
EKG system into an eigenvalue problem. The different 
solutions can be labeled by the value at the center of the 
first scalar field coefficient <f>\(r = 0) = (f>o, for which 
there is a set of eigenvalues of u>, 4>j>2{r = 0) = (fP Q , and 
Cj(r = 0) = Cq found to satisfy the boundary conditions. 
In practice, only the odd modes of the scalar field and 
the even modes of the metric functions are non trivial. 

As would have been guessed, the properties of oscilla- 
tons are similar to those of boson stars, and this has moti- 
vated the inclusion of oscillatons as models of galaxy dark 
matter |34j . However, they have not been as exhaustively 
studied as boson stars, and many of their properties may 
remain undiscovered. 



3. Stability 

Stability of boson stars has been studied both numer- 
ically and analytically 0, Restricting ourselves to 



the ground state configurations (also called 0-node solu- 
tions), it has been shown that there are stable equilib- 
rium configurations, generating the S-branch, which are 
stable against radial perturbations. That they are sta- 
ble is also indicated by the fact that S-configurations are 
indeed final states for a wide variety of initial configura- 
tions -including those in excited states- settle down onto 
at the end of numerical evolutions. Other initial config- 
urations either collapse to form a black hole or disperse 
away to infinity. 

On the other hand, oscillatons have been only barely 
studied and it is commonly believed that they are in- 
trinsically unstable. Though there is not an analytical 
proof for the intrinsic stability of oscillatons, numerical 
evolutions have shown that oscillatons can be classified 
into unstable and stable configurations, the latter be- 
ing called S-branch oscillatons, following the boson star 
nomenclature. S-oscillatons are stable against radial per- 
turbations, and they indeed play the role of final states of 
evolved systems made of real scalar fields 0. However, 
this is not a rigorous proof of the stability of oscillatons, 
and it is still possible that they have a long life time and 
decay very slowly [2^. 



C. The Newtonian limit 

Now, we are interested in the weak field limit version 
of the coupled Einstein-Klein-Gordon (EKG) equations, 
for which \\/ inG^ ^ (t, r)\ <C 1 and, correspondingly, 
\a 2 (t, r) - 1| < f and \a 2 (t,r) - 1| < 1. It turns out 
that the weak field limit of boson stars and oscillatons 
are quite similar. A common feature that simplifies the 
calculations is that weak equilibrium configurations have 
only one frequency, u> ~ m, whose corresponding Fourier 
mode is the dominant one in the case of oscillatons. 

As it was shown in 0, the weak field version of 
the EKG system is found through the standard post- 
Newtonian treatment, for both complex and real scalar 
fields. For this reason, we shall call them Newtonian 
boson stars and Newtonian oscillatons, respectively. 



1. Newtonian boson stars 

We start with the complex case following^. First, we 
express the scalar field and metric coefficients in terms 
of the Newtonian fields tp,U,A as 

^(t,x) = — L=e-* T ^(T,aO, (6a) 
V47rG 

a 2 {r,x) = 1 + 2U(t,x), (6b) 

a 2 {r,x) = 1 + 2A(t,x), (6c) 

where we have also introduced the dimcnsionlcss quan- 
tities r = mt, x = mr. Notice that U(t,x) and A(t, x) 
are real valued fields. 
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Next, we assume that the new fields arc of order 
0(e 2 ) <C 1, and that their derivatives obey the standard 
post-Newtonian rules d T ~ ed x ~ C(e 4 ). Therefore, con- 
sidering the leading order terms only, the EKG equations 
become 



id T ip 

d 2 x (xU) 
d x {xA) 



x 2 tpip* . 



(7) 

(8) 
(9) 



Notice that xp(r,x) is a typical complex Schrodingcr 
wave function, now coupled to U(t,x), which is the 
usual Newtonian gravitational potential. For this rea- 
son, Eqs. and (Et are named the Schrodinger-Newton 
(SN) system H |H[ll IH 13 ■ 



TABLE I: Quantities defined for the SN system. 



Density of particles 
Mass number 
Kinetic energy 
Gravitational energy 
Current of particles 



p(r,x) = tpi/)* , 
M(r,x)=J*py 2 dy, 
K{r,x) = -(1/2) /„" rd%(yip)ydy, 
W(r,x) = (1/2) /; pUy 2 dy, 
J(r,x) = (i/2)[ipd x iP* -ip*d x ip] . 



tions, in the weak field limit, reads 



2x 



[rdi{x^)-^di{xr)} , (12) 



d T A 2 = 2iA 2 - —i>d x ip . 



(13) 



2. Newtonian oscillatons 

The post-Newtonian approximation for oscillatons 
proceeds in a similar manner as for the complex scalar 
field. However, to take into account the oscillating na- 
ture of oscillatons we need to take extra fields, and then 
we should consider the expansions 

$M(r,s) = -^[e-^faaO+cx.] ' ( 10a ) 

a 2 (r,x) = 1 + 2U(t,x) + e- 2lT U 2 (T,x) +c.c.(lOb) 
a 2 (r,x) = l + 2^1(r,a;) + e" 2lT ^ 2 (T,x)+c.c.(10c) 

Notice that the fields ip(T,x), U(t,x) and A(t,x) have 
the same interpretation as in Eqs. JSJ). But, this time 
we have introduced two new complex fields U 2 (t, x) and 
A 2 (t, x) to take into account the oscillating nature of 
oscillatons. 

In any case, all fields are considered to obey the 
post-Newtonian rules depicted above. Therefore, the 
evolution equations for ip(r,x), U(r,x) and A(t, x) are 
Eqs. JJJ), (JHJ and In addition, an extra equation 

needed is 

d x U 2 = -x^ 2 . (11) 



3. Final remarks on the weak field limit 

Despite the presence of Eqs. @ and l|ll[l. it should be 
noticed that the dynamical evolution of both Newtonian 
boson stars and Newtonian oscillatons is dictated by the 
SN system only, Eqs. J7J and (JSJ. This only means that 
boson stars and oscillatons are quite similar at the weak 
field limit. 

So far, we have only used the {i, t} and {r, r} compo- 
nents of the Einstein equations to obtain Eqs. ©, © 
and jP) . The {t,r} component of the Einstein equa- 



Eq. i|12fl is just the conservation of probability as we 
know it from Quantum Mechanics. Eq. (|13f) only appears 
for Newtonian oscillatons, and represents the oscillatory 
behavior of the metric coefficient g rr . The rest of the 
Einstein components do not provide independent infor- 
mation. 

We should also mention here what we mean by the 
weakness of a Newtonian systems. The weakness of our 
configurations can be parametrized by the normalized 
value of the scalar field, and then we speak of the New- 
tonian limit if |^/4^G$( c • r )(^,^)| < 10- 3 [H. Any config- 
uration with a larger value should be treated using the 
fully relativistic EKG system. 



D. Properties of the Schrodinger-Newton system 

The properties of the scalar solitons, boson stars and 
oscillatons, can be more easily studied in the Newtonian 
limit, where we have a better understanding on the phys- 
ical processes involved. 

To begin with, we can unambiguously define many 
physical quantities that will help us to study the evo- 
lution and formation of weak scalar solitons. Some of 
them are listed in Tabic [I] These quantities are directly 
linked to their relativistic counterparts, whenever the lat- 
ter can be properly defined. For instance, the mass num- 
ber M(t, x) is related to the total mass Mr of the soliton 
configurations through 

f°° m 2 

M T = 4:iT p$(r, x)x 2 dx ~ —^M(t, x = oo) , (14) 
Jo m 

where p$ = — T* t in Eqs. (|TJ and J2J), for which we also 
find that p$ ~ (l/47r)m 2 m,p 1 ipifj* in the weak field limit. 
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1, Scaling properties 

Finally, we want to point out that Eqs. 1|7I8|) obey a 
scaling symmetry of the formpj 

{r, x, U, U 2 , A, i)} -> |A" 2 f, X~ l x, X 2 U, X 2 U 2 , A 2 i, A 2 ^} , 

(15) 

where A is an arbitrary parameter. 

This means that, once we have found a solution to 
the SN system, there is a complete family of solutions 
which are related one to each other just by a scaling 
transformation. Likewise, the quantities shown in Tabled 
obey the scaling transformation 

{p, M, K, W, J} -> |A 4 p, AM, X 3 K, X 3 W, A 5 j| . (16) 

Eq. (|15fl can be used to reduce significantly the space 
of possible equilibrium Newtonian configurations we need 
to study as follows. We will find solutions of the SN 
system for 'hat' quantities only, and then the 'non-hat' 
quantities will be obtained by applying the inverse of the 
scaling transformation 115fl . With this in mind, all quan- 
tities should be thought of as 'hat' quantities henceforth; 
we will write a 'hat' explicitly only when confusion can 
arise. 

Furthermore, we can always solve the 'hat' system tak- 
ing plainly that |?/>(0, 0)| = 1. Thus, the weak field limit 
condition is translated into a constraint on the scale pa- 
rameter as A 2 < 1CP 3 . 

Notice that the scaling transformation (|15|l does not 
apply for Eq. 113|) , but the latter indicates that A 2 (r, x) 
is, at least, A 2 -times smaller than the other Newtonian 
fields. In this respect, we can say that the g rr metric co- 
efficient is time-independent for both kind of Newtonian 
configurations. 



TABLE II: Eigenvalues of Newtonian equilibrium configura- 
tions. Shown also are their corresponding mass M, the 95% 
mass radius £95, and the kinetic K and gravitational W ener- 
gies. The precision of these numbers is in terms of the toler- 
ance we used in our shooting routine to solve the initial value 
problem, with the boundaries placed at least three times £95. 



Nodes 


7 


U(x = 0) 


M 


£95 


K 


W 





-0.69223 


-1.3418 


2.0622 


3.93 


0.47585 


-0.95169 


1 


-0.64793 


-1.5035 


4.5874 


8.04 


0.99037 


-1.9813 


2 


-0.63095 


-1.5811 


7.0969 


12.18 


1.4924 


-2.9850 


3 


-0.62081 


-1.6308 


9.5927 


16.35 


1.9850 


-3.9702 


4 


-0.61385 


-1.6670 


12.077 


20.54 


2.4706 


-4.9424 


5 


-0.60852 


-1.6954 


14.552 


24.74 


2.9520 


-5.9039 


12 


-0.58919 


-1.8060 


31.726 


54.40 


6.2333 


-12.463 



Although the system l|17fl has to be solved numerically, 
we can enunciate some analytical properties |l5t IT?! that 
simplify the numerical treatment. Eqs. 10 and (|17fl can 
be formally integrated up to 



4>{x) 

U(x) 
u 2 (x) 
A(x) 



1+ / y{l-y/x)(u-~/)<f>dy (18a) 
Jo 

M{x) 



U(0) 



yah dy 



-"2(0) - / y4> dy. 
Jo 

M(x) 



(18b) 
(18c) 
(18d) 



The relativistic condition of asymptotic flatness trans- 
lates into U(x — ► 00) = —M/x and u 2 (x — ► 00) = 0, 
from which we obtain 



E. Newtonian equilibrium configurations 

The SN system can be solved to obtain non-singular 
self- gravitating configurations. These so called equilib- 
rium configurations are of the form iI>(t, x) = e~ llT <p(x), 
for which Eqs. Q), ©, © and (fTTf) become the following 
system of ordinary differential equations 



(7(0) = -u 2 (0) 



(X(j>),xx 
(xTJ^j ^xx 
{U 2 ).x 

(xA), x 
- 2llT u 2 (x) 



2ar(Cr-7)0, 



X(j> , 

~X(j) 2 
2 j.2 



(17a) 
(17b) 
(17c) 
(17d) 



where U 2 (t,x) = e 

If we look for regular (0,^(0) = ^4(0) = 0) and finite 
configurations (<f)(x — > 00) — > 0), the SN system becomes 
an eigenvalue problem. For <fi(x = 0) = 1, there are 
unique values of 7, C/(0) and 1*2(0), for which the bound- 
ary conditions arc fulfilled. 



y<f> dy. 



(19) 



The numerical solutions are then found by using a shoot- 
ing method to adjust the value of 7, which is the only free 
eigenvalue, since the other two can be considered output 
values through Eq. (00. 

In Fig. n we show the equilibrium configurations for 
a 0-node and a 5-node configurations, and the required 
eigenvalues for different node configurations are shown in 
Table HH 

Some remarks are in turn. According to (|6a|) and (|10a|l . 
the fundamental angular frequency of the scalar field 
(j)(c,r) j g gj ven by tjj/ m = 1 -|_ ?y w ith 7 = A 2 7 Due 

to the weakness restriction, A 2 < 10~~ 3 , $( c - r ) has an an- 
gular frequency (in full units) of about co ~ m, but such 
that co/m < 1, as it is the case for gravitationally bound 
configurations [28| . 

From the Schrodinger equation Q, we find that the 
mass number (M), and the kinetic (K) and gravita- 
tional potential (W) energies are related through jM = 
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O-node equilibrium configuration 





<t>(x) 




U(x) 



















2 4 6 8 10 12 



x 

5-node equilibrium configuration 



2 




-1.5 

-2 I 1 1 1 1 1 1 

5 10 15 20 25 30 35 



x 

FIG. 1: Profiles of <j)(x), U(x) and U2(x) for a 0-node and 
a 5-node Newtonian equilibrium configurations, see text and 
TableETlfor details. 

K + 2W, while the total classical energy is Et = K + W. 
A common feature of any Newtonian equilibrium con- 
figuration is that, irrespectively of the number of nodes 
it may have, they all are virialized since iT/|VF| = 0.5, 
see Table HU Therefore, E T = (1/2) W = (1/3)jM < 0, 
which also indicates that all of the equilibrium configu- 
rations are gravitationally bound objects. 

F. Analysis of perturbations 

As we shall see later, our numerical simulations will 
be perturbed due to the discretization error, so we find 
illustrative to show how a first order perturbation model 
can predict extra modes that should appear in the evo- 
lution of equilibrium configurations. In this section, we 
restrict ourselves to the perturbations of 0-node Newto- 
nian equilibrium configurations. 

An important concept here is that of quasinormal 
mode. This is an oscillation mode which is characteris- 
tic of an equilibrium configuration which manifests when 
the system is perturbed. In relativistic studies, once the 
quasinormal frequency is known the mass of the equilib- 
rium configuration can be determined, thus it is possible 
to know the mass of the configuration a perturbed sys- 



tem is approaching to, without the need to follow the 
simulation until it settles down completely 0- We shall 
show that the same can be done for weak field solitons. 
First of all, we assume that 

tp = tp^+6ip \5ip\<l (20a) 
U = U {0) +SU \SU\<1 (20b) 

where -0'°' = <fi(x)e~ llT and U^°'(x) arc, respectively, 
the wave function and the gravitational potential calcu- 
lated for the unperturbed system (|1 7|> . Taking Sip = 
(p(x, t)e~ llt , the equations for the perturbations read 

id T cp = -l-d* x (xcp) + (tT(°) - 7) V + SUftla.) 
d 2 xx (x6U) = xW + tp), (21b) 

where we have neglected lower order terms. 

This system possesses similar scaling relations as those 
for the unperturbed SN: {SU,(p} {X 4 SU,X 4 ^}. Note 
that this indicates that the perturbations are suppressed 
faster than the unperturbed quantities as A — > 0. 

We then assume that the system <|21ll admits a station- 
ary solution of the form ip(x, t) = (pi(x)e~ la " r + ip 2 (x)e' l ' TT , 
where ip± and ip 2 are real functions and a is a real fre- 
quency. To first order, the resulting perturbations of the 
wave function, the energy density and the number of par- 
ticles read, respectively, 

8i>{x,T) = e"'^ T (<pie- iaT + y 2 e i<TT ) , (22a) 
6p(x, r) = <j> ( Vl + tp 2 ) (e" iCTT + e^ T ) , (22b) 

SM(t,x) = f Sp(y,r)y 2 dy, (22c) 
Jo 

and then we see that the perturbation of the gravitational 
potential admits the expansion SU(t, x) = 8u(x)(e~ laT + 
e lrTT ). Hence, the perturbation equations 12 II) are further 
reduced to 

dlxixipx) = 2x (V (0) -7 -crjtp! +2xSu(j) (23a) 
dl x (xcp 2 ) = 2x(u^ -7 + 0-) cp 2 + 2xSu^ (23b) 
d xx (xSu) = x(f) (ipi + tp 2 ) ■ (23c) 

The process followed to solve the system of equa- 
tions 112,'ill is the same as that used to solve the unper- 
turbed equations (|17f) . As before, we will restrict our- 
selves again to 'hat' quantities only in an implicit man- 
ner, so that we fix the central value of one of the functions 
arbitrarily to fi(0) = 1. The parameters ^2(0), Su(0) 
and a are free. Then, a shooting method is used to fix the 
values of the free quantities provided the boundary con- 
ditions SM(t,x — » 00) = ip\{x — > 00) = ip 2 (x — > 00) = 0, 
while at the same time we ask for regularity at the origin. 

The resulting perturbation profiles of ip±, (p 2 and du, 
are shown in Fig. [21 The angular cigcnfrcquency of the 
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perturbations is a — 0.2916, which coincides with other 
independent calculations ^1 1 . 





<PiM 




<P 2 M ■ 




5u x - 




SM(x) 





FIG. 2: Profiles of linear perturbations of the wave function 
denoted by (pi and <p2, and of the gravitational potential 8u 
for a 0-node equilibrium configuration, see Eqs. 121II . These 
perturbation modes correspond to a quasinormal frequency 
/ = 0.046. 



III. NUMERICAL TREATMENT 



In this section we give explicit details of the imple- 
mentation of a numerical code appropriate to solve the 
time-dependent SN system, Eqs. (0 and (jHJ. The issues 
covered are the discretization of the system of differential 
equations, boundary conditions and the accuracy of the 
solutions. 



The evolution 



The time-dependent Schrodinger equation is solved us- 
ing a fully implicit Crank-Nicholson scheme, which for 
equation J7J) reads 



1 + -HAt i/j 



/,n+l 



1 - -HAt i) n 



(24) 



The above solution means that the energy density 
and the gravitational potential of a 0-node configura- 
tion should oscillate with an angular frequency a when 
slightly perturbed (recall that the aforementioned quan- 
tities are strictly time-independent for an unperturbed 
configuration), and then / = a/(2n) ~ 0.046 would be 
the quasinormal frequency 2 of the Newtonian boson sys- 
tem when perturbed. It should be noticed that the quasi- 
normal frequency / is absent in the wave function per- 
turbation, but for this case, there are two perturbation 
frequencies, 7 ± a. 

The perturbation mode 5tp in Eq. I|20a|) is stationary 
since a is real. As it has been shown in^f|, this also 
means that the 0-node Newtonian configuration is stable 
against small radial perturbations. The stability analysis 
of excited equilibrium configurations will be performed 
with numerical simulations in the following sections. 



1 We should bo careful when compari ng values due to different 
scaling choices in the references. For ll5l , ta ke the first value of 
z>i in Table I, and then a = i>i/\/2. As for^lj: take the first value 
in Table II and use the conversion formula <r = o j V6 X 10 -2 . 

2 This should be compared with the values given by the for- 
mula (4.3) in 0, which in terms of 'hat' quantities reads 

-_ 7T 2 M 

When applied to a 0-node equilibrium configurations, it gives the 
value / = 0.033, which is at variance with our previous calcula- 
tions. Hereafter, we will refer to / = 0.046 as the correct value 
of the quasinormal mode of 0-node equilibrium configurations. 



being n the label of the time slice, At is the separation 
between time slices and the Hamiltonian H — — ^V 2 + £/ 
is written using second order centered finite differencing. 
This evolution scheme is found to be appropriate for the 
present problem since the approximation we use for the 
evolution operator is unitary, which allows one to pre- 
serve the number of particles M = J \ip\ 2 d?x [3fl| . 

The discretization of the Hamiltonian is second order 
accurate, and we handle the second term in the Laplacian 
as -j^- = 4-f%, being the last one a derivative with re- 

X ox ox 2 ' & 

spect to x 2 in order to avoid divergence at the origin. At 
the end, the finite differencing equation for the evolution 
of the wave function lf^|l is 

where p\ = Ar/4(Ax) 2 and P2 = Ar/(4Ax). The upper 
indices n still label time slices and lower indices label 
spatial grid points with j = 0,1. .A, where N labeling 
the last point of the grid. It is now defined a set of 
arrays that will help in the solution of this linear system 
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at each time step, 

aj = pi , (26a) 

Xj 

bj = -2pi--^+i, (26b) 
Ci = pi + — , (26c) 

Xj 

di = + i>U + (2 P1 + + i) ^ 

+ ("Pi - g) - (26d) 

for each grid point, with the special cases ao = 0, 
c N = 0, d = (V + ^ + i) Vo + (-Pi - fj) and 

4v = (-pi + f ) ^-i + (2pi + ^ +*) i>N- Then, 
the following linear system of equations arises 

6oV^o +1 +co^ +1 = do, 



aN-li>N + _ 2 + bN-l^N + -l + °N-li'N +1 = <^JV-1 , 

Such tridiagonal linear system is solved for the wave 
function ip n+1 at each time slice using a simple algorithm 
0. 



B. The Poisson equation 

In order to solve the Poisson equation we write it down 
in the form d 2 (xU) = x \ijj\ 2 and integrate inwards from 
the outer boundary. For this we use the Numerov algo- 
rithm, whose finite differencing expression for the gravi- 
tational potential reads 

(Ax) 2 

T njjn _ 2 n jjn _ n jjn , V I 

XjUj - zx j+1 u j+1 x j+1 u j+2 -t ^ X 

(^ +2 |^ +2 | 2 + 104 +1 |^ +1 | 2 + x?W\ 2 pt) 

This algorithm serves to integrate second order ordinary 
differential equations, and is locally sixth order accurate 
(see H3| for details). As we are solving the evolution 
equation (JJJ with a second order accurate differencing 
algorithm, the fact that we are solving the constraint (JSJ 
with a more accurate algorithm does not mean that the 
evolution itself should be better than second order. 

C. Boundary conditions 

Eq. (|27|l is to be integrated inwards and thus it is nec- 
essary to fix the value of the gravitational potential at 



the two outermost points of the numerical grid, i.e., we 
should deal now with boundary conditions. The easiest 
boundary condition for the gravitational potential is 

[/(zjv-i) = -M{x N - 1 )/x N ^ 1 , (28a) 
U(x N ) = -M(x N )/x N . (28b) 

However, this boundary condition has to be taken care- 
fully. Strictly speaking, jSEjl is only valid and consistent 
if the total mass is confined to the region x < xn-i (or 
cquivalently, M(xn-i) = M(xn)), since Eq. I)28|l is the 
form the gravitational potential takes in vacuum. 

In consequence, expression l|28f) is equivalent to impose 
the condition \ipij, xn)\ —* on the wave function. This 
is, for example, the shooting condition to find equilib- 
rium configurations we applied in Sec. Ill El Thus, we are 
forcing the system not to leave ever the computational 
domain, and then there cannot be outgoing waves at all: 
the outer boundary is a perfect reflective wall. 

The boundary condition is an important issue since 
there is no chance to apply the Sommerfeld boundary 
condition on the wave function ip as in the relativistic 
case 00 ■ because in the present situation the evolution 
equation is not hyperbolic. Moreover, if the boundary 
condition (|28|l were applied alone, it would only work 
for the equilibrium configurations -for which the num- 
ber of particles has to be preserved- but the evolution of 
other arbitrary systems, including those that eject mat- 
ter, would be incorrect. 

In order to avoid this, we implemented a sponge over 
the outermost points of the grid, which consists in adding 
an imaginary potential V(x) to the Schrodinger equation; 
the expression we use is 

i 

Vj = — Vo {2 + tanh [(xj — x c )/S] — tanh (x c /8)} , 

(29) 

which is a smooth version of a step function with am- 
plitude Vo and width 8. We choose this shape for the 
imaginary potential because in the absence of gravity, 
the wave function decays exponentially (see Sec. IIII C II 
below), and the smoothness is to diminish the effects of 
our finite differencing schemes. 

Notice that the definition of an imaginary potential 
makes physical sense and fits well into the SN system. 
This is seen if we recalculate the conservation of proba- 
bility which now reads (see Eq. dJ) 

d T p + y x WW) ~ rd 2 x (xyj)] = -2iVp . (30) 

The source term on the r.h.s. is always proportional to 
the density of particles, and the minus sign warranties 
the decay of the number of particles at the outer parts of 
our integration domain, that is, the imaginary potential 
behaves as a sink of particles. 

One last remark. Our original SN system is not 
physically related to the sponge region. Then, we ad- 
just the parameters of the sponge so that it only cov- 
ers some points at the outer part of our grid; actu- 
ally, we have fixed the spatial range of the sponge to 
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2(x N 



The values of A s 



and x 7v 



are fixed by hand and then we calculate x c and set 
S = A spongo /fO. Notice too that Vn = —iVq. In this 
manner, we shall call physical region to the grid points 
that accomplish Xj < (xjv — A spongo ), and this will be 
our region of physical interest. 



1. Toy model: Imaginary 1-D square well potential 

However, the sponge <|2"§1) is not a perfect absorber, and 
to see how it works in reality we consider the following toy 
model (a more general discussion can be found in 32] ). 

We simplify the system and consider that at the last 
points of our numerical grid, where the sponge is im- 
plemented, the gravitational potential can be neglected. 
Then, our physical system at the last points would be 
similar to the case in which a Schrodingcr wave function 
is scattered back off by an imaginary square well poten- 
tial of the form V(Q < x < x c ) = — iVo and V(x) = 
elsewhere. 

As in usual quantum mechanics examples, the radial 
Schrodinger wave function xip(r, x) = u r (x)e~ lET has the 
form 



Re ~ikx 



u r (x) 



p e ik'x + Q e 





-ik x 



x < , 
, < x, < I 
£ < x , 



(31) 



where k = \J~EJVq and k' 2 = i + k 2 are the wave numbers 
outside and inside the well potential, respectively. The 
presence of imaginary number i reveals the imaginary 
nature of the well potential. 

The coefficients R, P, and Q are already normalized 
with respect to the incident wave e tkx , while the spatial 
coordinate is normalized in the form x = x\/2Vq and 
then £ = A sp0 ngeV2Vb. Notice that we have added the 
condition of perfect reflection at the right boundary x, = 
£, as it is the case for our numerical system. 

Using the continuity and differentiability conditions for 
the wave function at the boundaries x = and x = £, we 
find the following expression for the reflection coefficient 



\R\ 



iksin(k'£) + k'cos(k'L) 



i(k'£) - k' cos(fc'L) 



(32) 



This formula is also valid for a real well potential under 
the same boundary conditions, and for such case it gives 
the obvious solution |i?| 2 = 1. 

A graph of the reflection coefficient |i?| 2 as a function 
of (k, £) is shown in Fig. [31 where we see that |i?| 2 is highly 
suppressed for k 3> 1 and I 3> 1. The former condition 
means that more modes can get into the sponge, while 
the latter one means that the sponge region is sufficiently 
large as to let the modes die out inside the sponge. 

Going back to our original system, it is then obvious 
that low energy modes will always be reflected in some 



|R|»2 




FIG. 3: The reflection coefficient \R\ 3 , Eq. , as a func- 
tion of k and I, for the linear toy model shown in text. The 
reflected number of particles is suppressed for large values of 
both parameters k, I. 



amount and they will always contaminate what we called 
the physical region. However, the reflected amount of 
matter is not as large as in the toy model above, since 
the imaginary potential (|29|l is a smooth function. Then, 
the depth and the size of the imaginary potential can be 
adjusted to let most of the modes enter the sponge. 

Fortunately, because of the fitness of the imaginary 
potential into the SN system, the sponge region could 
be larger than the region of physical interest without re- 
sulting in an erroneous evolution of the SN system. The 
appropriateness of the sponge is analyzed in Sec. IIV Al 



D. Accuracy 

In the relativistic case, the EKG system is redundant 
and one of the Einstein equations can be used to mea- 
sure how accurate the numerical solution is. For in- 
stance, in|{J, the {t, t} and {r, r} components of the Ein- 
stein equations were solved and then it was determined 
whether the {t, r} component was accurately satisfied. 

Following the relativistic case, Eq. I|3U|I gives a criterion 
to measure the accuracy of our numerical code for the 
complete grid of integration, since it is an equation our 
system should accomplish independently of the potential 
involved. As said before, such equation is the weak-field 
constraint equation of the relativistic system. Then, we 
define our accuracy parameter to be: 

[3 := d T p + ^- [V^(:#*) " 4>*dl{x4>)} + ZiVp . (33) 

In the continuum limit, (3 = holds for an exact solution, 
thus the magnitude of this quantity will tell us how close 
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our numerical solution is to the exact one. 

However, it should be noticed that our accuracy pa- 
rameter obeys the scaling transformation (3 = A 6 /3. So, 
the value of /3 should be interpreted carefully otherwise 
we could claim that the most accurate physical solution 
is that with the smallest value of A, which cannot be true 
(accuracy cannot be just a matter of scaling). 

In most cases (see for instance Sec. IIV D|) . it would 
be more useful to define a A-independent parameter. A 
simple and straightforward definition we can think of is 
a relative error parameter of the form 



A/3 



1 1 St A 

EillAI 



(34) 



Here, Pi is each of the terms that appear in the momen- 
tum constraint, whether relativistic or non-relativistic. 
For instance, pi represents each of the terms we need to 
calculate our accuracy parameter P in Eq. (|33|l : seeQ for 
the definition of a P for a relativistic system. Clearly, 
A/3 = A/3 for the SN system, and we will indicate its 
value whenever the interpretation of P is confusing. 

As a final note, we should mention that the accuracy of 
the evolutions depends on the ratio At/(Ax) 2 , the lat- 
ter should be less than unity in order to have a reliable 
evolution of the Schrodinger equation, see|3(]]. However, 
it should be noticed here that, according to the scale 
transformation (fT5|> . At/ (Ax) 2 = Ai/(A£) 2 , and then 
the computational effort is the same for the original con- 
figuration as for the properly sized one. Actually, the 
condition Ar/(Ax) 2 < 1 is a very restrictive one, and 
the responsible for the very large runs we need in order 
to get a reliable evolution at late times. 



IV. RESULTS 

In this section, we show representative runs of the dif- 
ferent issues discussed in the previous sections. The nu- 
merical method to evolve the SN system is tested thor- 
oughly about the effect of the boundary conditions, accu- 
racy and convergence of the numerical results. In brief, 
the tests show that the numerical results arc reliable even 
in long runs, as far as the boundary conditions are set up 
appropriately. 

Another test considered is the evolution of 0-nodc 
equilibrium configurations and their perturbations. The 
numerical evolution reproduces the semianalytical results 
of the previous sections and shows that 0-node equilib- 
rium configurations are intrinsically stable. On the other 
hand, it is shown that the numerical method preserves 
the scaling properties of the SN system even in the cases 
of non-equilibrium solutions. Also and for completeness, 
the equivalence of the (relativistic) EKG and the SN sys- 
tems is studied in the weak field limit and found to be 
correct within the numerical accuracy. Finally, excited 
equilibrium configurations are found to be intrinsically 
unstable and that they all decay to 0-node solutions. 



A. Boundary effects 

Being the boundary conditions an issue of particular 
importance in this work, it is important to determine the 
reliability of the sponge for the purposes of this paper, 
since part of the initial mass could be ejected to infinity 
and forced to interact with the imaginary potential at 
the outer boundary. 

The form in which we test our boundary condition is 
as follows. The ideal numerical evolution is that in which 
the numerical boundary is very far away from the phys- 
ical boundary, x p <C xn, so that no scalar matter has 
reached xn for the time range the numerical evolution is 
followed. In this manner, we can assure that the solution 
inside the physical region < x < x p is the correct one. 

On the other hand, let us suppose that x v ~ xn- and 
that a sponge is implemented on the region x p < x < 
xn- If the numerical evolution in the range < x < x p 
coincides with that of the ideal case above, then we say 
that an appropriate sponge was implemented. 

With this in mind, we have numerically evolved the ini- 
tial Gaussian profile shown in Fig. 0] which is of the form 
■0(0,2;) = exp[— (x/2) 2 ]; also shown is the corresponding 
metric function A(0,x), see Eq. (|10c|) . We have marked 
the radii £95 (the so called 95% mass radius) and x max , 
where the latter indicates the radius at which A(t, x) 
reaches its maximum value. In the relativistic systems, 
Xmax has been a good quantity to follow the evolution of 
scalar systems, see [a. l27j. 

The time and space resolutions of the numerical evo- 
lution were, respectively, At = 3 x 10 -3 and Ax = 0.08. 
Also, the numerical and physical boundaries were xn = 
960 and x p = 40, respectively, and no sponge was imple- 
mented. We shall call this case Run I. 




FIG. 4: Initial profiles of the wave function 0(0, a;) = 
cxp [-(x/2) 2 ] (left axis) and the metric function A(0, x) = 
M(0, x)/x (right axis) of Run I. x max indicates the position 
of the maximum value of A(0,x), and £95 is the 95% mass 
radius. 

In Fig. [5] we show the evolution of a; max and £95 for 
Run I. We notice that x max oscillates and approaches to 
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a finite fixed value, while X95 grows almost linearly for 
all the time the evolution was followed. 

Roughly speaking, x max gives us an estimate of the 
radius at which most of the gravitational interaction is 
contained in, while X95 can be seen as a tracer of the 
ejected mass that escapes to infinity. Thus, we can say 
that a finite sized configuration is formed, and that the 
ejected scalar matter is far from reaching the numerical 
boundary. That is, the solid curves in Fig. [S] and arc 
what we would see in the ideal numerical evolution. 

Shown in Fig.[S]is the mass contained inside the phys- 
ical region < x < 40 as the evolution proceeds. For 
the case of Run I (solid curves), we notice that part of 
the mass leaves the physical region (r ~ 300) but a little 
bit of it returns (r ~ 500) and leaves again (r ~ 600). 
This return of scalar matter cannot be attributed to re- 
flection at xn since no matter has reached the numerical 
boundary. But, this effect is due to matter that leaves 
the physical region with a velocity less than the escape 
velocity of the system. 

Being Run I an example free of noise from the nu- 
merical boundary, we will compare it with other runs in 
order to test the reliability of a sponge as a boundary 
condition. For sake of simplicity, the width of the imagi- 
nary potential was fixed to 5 = (xn — x p )/10 in the runs 
described next. 




100 200 300 400 500 600 

X 

FIG. 5: Evolution of the two radii x max (left axis) and X95 
(right axis) for Runs I, II, III and IV; see text for details. 

Run II has the same time and space resolutions as Run 
I, but the values of the physical and numerical boundaries 
are, respectively, x p = 40 and xn = 240, and the sponge 
depth (see Eq. J25JI) is Vq = 55.0. For this case, X95 
reaches a maximum value but never reaches the value 
of the numerical boundary; this means that the ejected 
matter leaves the physical boundary (from r ~ 100 on) 
and is absorbed by the sponge. 

On the other hand, the value of x max (see Fig. [SJ) again 
oscillates and approaches to a finite value, but we note a 
shift in the oscillations with respect to Run I. The reason 
for this can be found if we look at the plot of the mass 
number M(t, x = 40) in Fig. 
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FIG. 6: Mass number M(r,x = 40) (left axis) for all Runs I 
to IV. Also shown is the relative mass difference (right axis) 
of Runs III and IV with respect to Run I. For Run IV, the 
mass deviation is almost unnoticeable; see text for details. 



First of all, the reflection of matter is noticeable around 
t ~ 200, which is due to the sharpness of the sponge, 
recall that Vq = 55.0. Also, the size of the physical region 
is very small, since the sponge absorbs matter that should 
have returned to the physical region because it did not 
have the required velocity to escape to infinity, as it is 
shown by Run I. At the end of the run, there is a mass 
difference of about 5% between Run I and Run II in 
Fig. El which is sufficient to make Run II evolve apart 
from Run I at late times. 

Run III uses the same parameters as Run II, but now 
Vq = 1.0. In general, the results are qualitatively the 
same as for Run II, but the amounts of reflected and 
absorbed matter are drastically reduced now that the 
sponge is smoother. As seen in the corresponding X95 in 
Fig. the ejected matter can travel a longer distance 
inside the sponge region and some of it can return to 
the physical region without being completely annihilated. 
The mass difference with respect to Run I at the end of 
the run is around 0.1% (see Fig.|H}, and then the oscilla- 
tions of x max arc not shifted noticeably with respect to 
those of Run I. 

Run IV has the same parameters as Run III, except for 
the physical and numerical boundaries now at x p = 80 
and at x n = 280, respectively. That is, the sponge region 
has the same size in Runs II, III and IV. The latter is the 
most similar to Run I, being the mass deviation in Fig. El 
less than 0.001% at the end of the simulation. This is 
because the physical region is larger than in Run III, and 
then more ejected matter can return to the region x < 40 
without being affected by the sponge. 

All the results presented here show that with an ap- 
propriate sponge we can obtain the same results as if the 
numerical boundary were very far away, with the advan- 
tage that Run IV needs less numerical efforts than Run 
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I 3 . 

Finally, we show in Fig. a comparison of the final 
profiles of x p corresponding to Run I and Run IV. As 
it should be for a good boundary condition, the profiles 
are quite similar inside the physical region, which is our 
region of interest. Moreover, it is seen that the expected 
behavior |^>(r, xn)\ — *■ is achieved, which is the bound- 
ary condition compatible with Eq. (|28(l . 




50 100 150 200 250 300 350 

X 

FIG. 7: Comparison of the final profiles of x 2 p(r = 600, x) 
corresponding to Run I and Run IV, where we have marked 
the physical and sponge regions for the latter. The profiles 
are quite similar inside the physical region, but the profile of 
Run IV rapidly vanishes inside the sponge region, and then 
we are accomplishing the boundary condition \ij)(T, xn)\ — * 0. 

In conclusion, the implementation of a good sponge 
should follow, in general, the indications found in 
Sec. IIII C II a smooth and large sponge will make it well. 
However, the examples presented in this section gives 
more details about how smooth and how large a sponge 
could be in order to have a reliable run and a low numer- 
ical effort. Also, we learned that a good sponge should 
be accompanied by a sufficiently large physical region. 



B. Accuracy and Convergence 

Once we have shown the reliability of the boundary 
conditions, we can perform accuracy and convergence 
tests. 

To begin with, we show in Fig. [S] how ||/3||2 evolves in 
time; its value is less than 10~ 8 for Runs I, 77, III and 
IV. The values of ||/3||2 are the same because the ratio 
At/(Ax) 2 is also the same for all runs. It is clear here 
that all deviations from Run I are due to matter reflected 
by the sponge. 



3 The runs were followed up to f = 600, but recall that the physical 
dimcnsionless time is such that r > 10 3 f ; thus, as far as we know, 
the runs in this section are the largest reported in the literature. 
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FIG. 8: (Top) The value of |fS|] 2 for Runs I, II, III and IV. 
In all these cases, A/3 < 10 -6 . (Bottom) Evolution of <r max 
(left vertical axes) for the same case as in Run IV, but with 
a fixed time step At = 2 x 10 -4 and four different spatial 
resolutions Ax = 0.02, 0.04, 0.08 and 0.16. It can be seen 
that the numerical solution converges as Ax — > in phase and 
amplitude. Also shown (right vertical axes) is the difference 
in Zmax between runs with spatial resolutions Ax and 2 Ax; 
the difference is four times smaller if the spatial resolution is 
doubled, see text below. 



Another test we consider important is that of conver- 
gence, that is, whether the solutions of our numerical 
scheme approach the exact solution as we increase the 
spatial resolution. To investigate this, we perform four 
different runs for the same system as in Run IV with a 
fixed time step At = 2 x 10 -4 , and four spatial resolu- 
tions Ax = 0.02, 0.04, 0.08 and 0.16; the corresponding 
evolutions of £ max are shown in Fig. [5] 

Only the coarsest run (Aa; = 0.16) deviates to second 
order in phase from the other three, and it was found that 
such coarse resolution is out of the convergence regime af- 
ter certain finite time of evolution; such deviation should 
be attributed to the discretization used and not to re- 
flected matter. Moreover, it can be noticed that the so- 
lutions indeed converge to the solid line as the spatial 
resolution is increased. 

We also show the deviations when we compare the runs 
by pairs, f^ 2Ax ^ - f( Ax \ where is any function of 
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our system that was solved using a fixed spatial resolu- 
tion A. For the case of £ max , we notice in Fig. |SJthat 
the deviations are four times smaller when the spatial 
resolution is doubled. 



C. Equilibrium configurations in the ground state 

Another important test is provided by the existence 
of stationary solutions of the SN system of equations 
which were found in Section III El In short, these so- 
lutions are obtained by assuming that the time depen- 
dence of the Schrodinger field is of the form ip = e 47T </>(r), 
which in turn implies that both the gravitational poten- 
tial U(x) and the probability density p = \ip\ 2 are time- 
independent. 

In Fig. El h is shown the evolution of Ke(ip(r, 0)) for a 
0-nodc equilibrium configuration. The boundaries were 
located at x = 50 with a time resolution of At — 0.001 
and a spatial resolution of Ax = 0.02. It is observed that 
the wave function oscillates harmonically, as its Fourier 
transform (FT) shows a unique harmonic mode with an 
angular frequency 7 = 0.697, in good agreement with the 
eigenvalue solution found in Section Hi El 

As was discussed in Section IIII CI the boundary 
condition (|28[) imposed on the gravitational potential 
makes the energy density vanish at the outer boundary. 
As there is not any violent collapse or explosion in 
this system, such boundary condition is appropriate for 
0-node equilibrium configurations. Fig. [5] also shows 
that the numerical code is stable and, within numerical 
precision, reproduces the expected analytical results. 



1. Observing extra modes 

The evolution of the 0-node configuration was ob- 
served carefully in order to look for the modes predicted 
by the perturbation theory described in Sec. Ill Fl For 
this, we searched for the perturbations of the energy 
density as they can be easily isolated, see Eqs. (|2"2|) . 
We show in Fig. ^] the oscillations on the energy den- 
sity of the system due to the perturbations introduced 
by the discretization of our grid. Actually, we are 
comparing three runs with three different resolutions 
Ax = 0.04, 0.08, 0.16 with the time resolution fixed at 
At = 10" 3 . 

It is noticed that the perturbations with the smallest 
amplitude corresponds to the run with the finest spatial 
resolution, as it can be seen from the FT of the oscil- 
lations also shown in Fig. This plot also indicates 
that the oscillation frequency / ~ 0.046 is the same de- 
spite the resolution used, which coincides with the value 
found by solving the perturbation equations for the 0- 
node configuration in Section Til Fl Last but not least, 
the amplitude of the perturbations also show that the 
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FIG. 9: (Top) The evolution of Re(tp(r,0)) for a 0-node 
equilibrium configuration. (Bottom) The DFT of the evolved 
solution shows a unique harmonic mode for an angular fre- 
quency of 7 = 0.697, which coincides well with the eigenvalue 
problem of Section Hi El 

numerical code is second order convergent: for a fixed 
At, the amplitude of the perturbations are four times 
smaller if we double the spatial resolution, in concord 
with the results of Sec. IIVBI 

For further comparison, we show in Fig. llll the profiles 
of the perturbations as obtained from the evolved sys- 
tems and the perturbation equations H23|) . Although the 
correspondence is not exact, the results still suggest that 
the perturbations in Figs . l2l and 1 1 01 are of the same kind. 

A general conclusion of this section is that, as antic- 
ipated from the perturbations analysis, 0-node equilib- 
rium configurations are intrinsically stable under small 
radial perturbations. 

D. Scaling relations 

For the scale invariance (|T5|l to be really useful, we have 
to be sure that the initial scaling relation is preserved by 
the numerical code during the evolution of the SN system. 
Also, we would like to test the post-Newtonian approxi- 
mation made for the relativistic system, i.e., whether the 
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FIG. 10: (Left) The central value p(r, 0) at different resolu- 
tions Ax = 0.04, 0.08, 0.16 using a fixed value of At = 10~ :! . 
The plots show that the numerical code is second order con- 
vergent. It can be observed that for the finest resolution there 
is an initial noise due to the fact that the quantity At /{Ax) 2 is 
bigger than in the other two cases. (Right) The FT of p{j, 0) 
shows the first peak which we associate to the quasinormal 
frequency at / ~ 0.046, see Sec. Ill Fl 



SN system is the weak field limit of the EKG equations. 
That the latter is true has been shown before in the case 
of boson starsj^j, but we want to show it also for the 
case of oscillatons. 

To test both features above in one stroke, we show 
the numerical evolution of a 0-nodc Newtonian oscilla- 
ton which was perturbed so that its mass was increased 
roughly by 20%; such initial profile was used to feed 
the relativistic numerical scheme used in@, following 

Eqs. (frn|>. 

The scale invariant quantities we used for the SN sys- 
tem were Ax = 0.04, f = 3.0 x 10~ 4 , with the physical 
and numerical boundaries at x p = 80 and at xn = 160, 
respectively. The numerical evolution of the SN system 
was followed up to a time f = 27. 

On the other hand, to set up the initial profile for the 
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FIG. 11: The normalized perturbation profiles of the energy 
density (left axis) and the gravitational potential (right axis) 
of a 0-node equilibrium configurations. Shown here is a com- 
parison between the evolved profile shown in Figs. [5] and 1101 
and the solutions obtained in Section fll Fl Although there is 
not a complete correspondence, the similarity of the graphs 
suggests that we indeed see the same kind of perturbation in 
both cases. The spatial resolution of the evolved profile was 
larger than shown in the plot. 



relativistic case, the EKG system, we take the scaling pa- 
rameter A = 0.032. Because of the scale transformation, 
the spatial and time resolutions used were Ax = 1.25 and 
At = 0.0125, respectively, while the physical boundary 
was fixed at x p = 2500. The relativistic run was followed 
up to t = 25000, that is, just 10 crossing times. 

The resulting graphs of the energy density p(r, 0) for 
both the (post-Newtonian) SN and the (relativistic) 
EKG systems are shown in Fig. To begin with, a 
simple comparison of the scaled and non-scaled quanti- 
ties tells us that the SN evolution is simpler and needs 
less computational effort when evolving weak field con- 
figurations. Also, we see from the plots that both evo- 
lutions are quite similar and then we can be confident 
that the SN system is also the weak-field approximation 
of oscillatons. 

The discrepancies seen in the graphs, which are less 
than 4%, can be attributed to the numerical error of 
the relativistic system (the system is too weak and the 
metric functions are pretty close to the values of the 
flat space, which makes them difficult to evolve accu- 
rately), and to the fact that we are at the edge of 
validity of the weak-field approximation. As it was 
first discussed in[2^|, Newtonian oscillatons are valid if 
V87rG|$| = 2A 2 < V2 x 10~ 3 , and the scale parameter 
used in our example is just above this limit. Had we 
taken a weaker configuration, the EKG solution would 
have been even less accurate. 

In order to support the latter statement, we show 
in Fig. Ua the relative error A/3, Eq. (|34|l . for the 
relativisticlal and non-relativistic runs. Thus, we con- 
is more accurate 



elude that the SN result 
than the EKG one (A/3 - 



(A/3~ 
; 10~ 3 ). 
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FIG. 12: (Top) The numerical evolution of the maximum of 
the matter density p m ax{T~) using the relativistic EKG system 
and a properly sized SN system. The initial scalar profile was 
a perturbed Newtonian oscillaton with 20% of mass excess, 
see text for details. (Bottom) Comparison of the relative vi- 
olation of the momentum constraint A/3 for both the EKG0 
and the SN system, see Eq. 1341 . It can be seen that the 
Newtonian run is more accurate than the relativistic one; the 
reason for this difference is that the evolved system is too 
weak for the relativistic code. In both figures, the axes used 
for the relativistic system are labeled (EKG), while the label 
(SN) appears for the non-relativistic system. The output data 
files for both systems are plotted without scaling, instead the 
scale transformation relating the runs can be calculated from 
the ranges shown in the axes. A = 0.032 for the case shown 
in here. 



We have also evolved the SN system for different initial 
configurations related initially by a scaling transforma- 
tion. An example of two Gaussians related initially by 
ipo — 4-00 is shown in Fig. 1131 where we plot the resulting 
gravitational potential. It can be seen that the scaling 
transformation is preserved always. In fact, for all cases 
we studied, the scaling transformation (|15fl was obeyed 
very accurately all along the evolution 

The scaling transformation of the perturbations in the 
energy density can also be seen in the evolved systems. 
In Fig. ll4l we show again the perturbations of the configu- 
ration that appears in Figs. 151 and II II (left axis), and the 
corresponding ones for the scaled system with A = 0.1 
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FIG. 13: The minimum of the gravitational potential for the 
initial configuration %[>(x,6) = e"^ 2 ' 75 ' is shown. Super- 
posed it is also the rescaled evolution for the corresponding 
configuration using A = 2, that is ip(x,0) = 4.0e" (:c/1 ' 375) " . 
The data files were not modified, and the ranges shown in 
the axis confirm the scaling properties of the SN system. For 
this case, |Z7 m i n — A 2 ?7 m m| < 2 x 10~ 6 , i.e., the difference is 
unnoticeable in the plot. 

(right axis). As stated in Sec. Ill Fl it is easily seen that 
the perturbations are related by the scaling transforma- 
tion Sp = X 4 Sp. 
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FIG. 14: Scaling transformation of the density perturbations 
Sp(r, 0) shown in Figs. |5] and 1111 The scale parameter used 
was A = 0.1. The data files were not modified, and the ranges 
shown in the axis confirm the scaling properties of the SN 
system. It is seen that the relation for the perturbation Sp = 
\ 4 Sp is satisfied, see Sec. Ill Fl since \5p — A 4 <5/3| < 4 x 10~ 9 . 



E. Equilibrium configurations in excited states 

Excited equilibrium configurations (also called n-nodc 
configurations) are also stationary solutions of the SN 
system, and the aim of this section is to investigate 
whether they arc stable. What we find here is that, in 
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general, all excited configurations are intrinsically unsta- 
ble and decay onto 0-node configurations by emitting 
scalar matter. 

As a representative example of such decay, we show 
the numerical evolution of a 1-node equilibrium configu- 
ration in Fig. El Even though this excited configuration 
is initially virialized (if/|W^| = 0.5) and only perturbed 
by means of the discretization error, it ejects mass and 
settles down onto a 0-nodc equilibrium configuration, as 
it can be seen from the plots of the profiles of x 2 p as the 
evolution proceeds. 

At later times, the energy density of the evolved 1- 
node system oscillates with frequency / = 0.0976, which 
means that it will settle down onto a 0-node equilibrium 

configuration with a scale parameter A = \J f / f = 1.457, 

where / is the quasinormal mode found in Sec. Ill Fl That 
is, the final state is that with M = 3.0 and x max = 1.77. 
This last fact can be seen in a plot of M vs x max , as in 
Fig. El where the migration path of a 1-node system 
can be followed until it ends up at a 0-node equilibrium 
configuration. 




FIG. 15: (Top) The virialization K/\W\ of a numerically 
evolved 1-node equilibrium configuration. The system is ini- 
tially virialized but it is nonetheless unstable. (Bottom) Some 
profiles of x 2 p for a 1-node equilibrium configuration at dif- 
ferent times of its evolution. At the end of the run, the system 
has lost its node and evolves towards a 0-node configuration 
represented by the thick solid curve, see also Fig. ED 



The migration paths of different excited equilibrium 
configurations are also shown in Fig. 1161 Without any 
added perturbation (apart from the discretization error), 
they all decay and migrate to a 0-node solution. The 
latter are represented by the solid curve drawn by the 
formula 



M = M 



0-node / X 



„0— node 



(35) 



with M°- nodc 



2.0622 and x^™ dc = 2.58. This last 
formula appears from the fact that the scale parameter 
is given by A = M/M = x max /x max , sec Eq. (JTSJ). 




Migration 



5-node ■ 

4-node ■ 
3-node 

2-node ■ 

1-node ■ 

0-node ■ 



25 



FIG. 16: Plots of the evolution of the mass number M and 
2; ma x for different (unperturbed) excited equilibrium config- 
urations. The solid curve marks all possible 0-node systems 
related through a scaling transformation, that is, those rep- 
resented by the formula (1351 . 



V. GRAVITATIONAL COOLING 

In this section, we investigate an interesting is- 
sue which arises in the formation of gravitationally 
bounded objects, in particular, the scalar objects of this 
manuscript. 

It was discovered numerically 3 that, for the SN sys- 
tem, an arbitrary initial configuration settles down into a 
stable configuration through the emission of scalar mat- 
ter, with the property that the ratio K/\W\ approaches 
and oscillates around 0.5 at late times, which means that 
the system is around virialization. 

This so-called gravitational cooling is so efficient, that 
allows the virialization of overmanned systems for which 
K/\W\ > 1 initially. This is to be compared to the in- 
ability of the violent relaxation process to virialize such 
systems. 

For the rest of this section, we study with more de- 
tail the gravitational cooling of arbitrary scalar systems 
whose evolution is dictated by the SN system. As we 
shall see, the results confirm the ability of the gravi- 
tational cooling to virialize all possible initial configu- 
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rations, even if they are too overwarmed. This section 
complements the relativistic studies in[(|. 

To begin with, we evolve the same configuration pre- 
sented in 0|, but properly scaled in the form 



tp{0,x) = 4.5e 



-(x/2.1) 



[1 + 0.5 sin(7rx/0.15)] (1 + i) ;(36) 



6 x 10 4 and x p = 1500 as 



that is, we are using the scale parameter A = 0.01. The 
spatial and time resolutions were Ax = 5 x 10 -3 and 
Af = 1.25 x 10 -5 , respectively; and the run was followed 
up to f = 6 while the physical boundary was set at x p = 
15 (which correspond to r 
ini). 

In Fig. El we show the same quantities shown in Fig. 3 
which are the virialization K/\W\ and the total en- 
ergy of the system E = K + W; the corresponding graph 
of the mass is shown in Fig. 1 181 and will be discussed later. 
Even though the results agree qualitatively, we find that 
the values of the mass and the total energy at the end of 
the run arc not the same. 
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FIG. 17: The virialization K/\W\ (left axis) and total energy 
E = K + W (right axis) for the initial distribution (EHJ. The 
graphs are scaled for an easy comparison with Fig. 3 infl- 



The system seems to settle down onto an equilib- 
rium configuration with a (non-scaled) mass of around 
M = 0.45. But, this value is within the relativistic 
realm, since for Newtonian equilibrium configurations 
M < 0.08. Thus, we consider it would be more appropri- 
ate to evolve the system <|36[) using the relativistic EKG 
equations. 

We made a numerical run of the EKG system for a 
real scalar field, using the following initial profiles for the 
scalar field and its time derivative, respectively, 



/^$ (r) (0,x) = 2A 2 Re(i/>(0,x)) ; 
'kq— — (0, x) 



dr 



2A 2 Inh>(0,a;)), 



(37a) 
(37b) 



as follows from Eq. l|10a|) and the post-Newtonian rules. 
The spatial and time resolutions of this relativistic evo- 
lution were Ax = 0.5 and At = 0.01, respectively, while 
the numerical boundary was set at xn = 1500. 



For comparison, we plot in Fig. ^| the total mass of 
the SN system with that of the EKG system. Initially, 
both systems have the same mass, but the path followed 
by each system is different as the evolution proceeds. As 
we mentioned before, we believe that the true final state 
is that of the relativistic system, which in this cases cor- 
responds to a relativistic oscillaton. 

At this point, we would like to mention that the com- 
parison between violent relaxation and the gravitational 
cooling is not appropriate in this example. As Fig. [TBI 
shows, the initial configuration is so massive that it can- 
not disperse away because of the intervention of strong 
gravitational effects fully accounted by the Einstein equa- 
tions, despite the large initial ratio K/\W\ = 1.4. Actu- 
ally, the inclusion of relativistic effects can prevent the 
appearance of the gravitational cooling and make the sys- 
tems collapse into a black hole|i|. 
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FIG. 18: Comparison of the evolution of the total mass for 
the same Newtonian system shown in Fig. 1171 (see Eq. 1361 '). 
and the relativistic system represented by Eqs. 1371 . see text 
for details. As in Figs. 1121 and 1131 the output data files for 
both systems are plotted without modification, and the scale 
transformation relating the runs can be calculated from the 
ranges shown in the axes. A = 0.032 for the case shown in 
here. 

Nevertheless, we have indeed found that the gravita- 
tional cooling is a very efficient mechanism for Newtonian 
systems (for the relativistic case see 0,0) even in the case 
in which K/\W\>1. 

To investigate this, we have made runs with an initial 
Gaussian profile of the form i/j(0,x) = ipoe - ^/ 2 ^ tak- 
ing different values values of i/V In Fig. ^3 we show 
the graphs of the virialization coefficient K/\W\ and the 
corresponding path in a M vs x max plot of the evolved 
Gaussian profiles. 

For the value ip = 2.0 (K/\W\ T=0 = 0.265), the sys- 
tem rapidly virializes, and, from r = 1500 afterwards, 
it has settled down onto a 0-node equilibrium configura- 
tion corresponding to a scale parameter A = 2.33. For 
the values Vo = 1-0 (K/\W\ T=0 = 1.06) and ip = 0.98 
(if/|W| T= o = 1.104), it is noticed that the systems are 
approaching to 0-nodc equilibrium configurations, but it 
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FIG. 19: Evolution of different overwarmed and overcooled 
Gaussian profiles of the form ip(0, x) = i/joe"^ 2 ' . Shown are 
the cases ip = 2.0, 1.0,0.98,0.93,0.92. (Top) Independently 
of the initial value of if/|W|, the scalar systems always col- 
lapse back to an equilibrium configuration. The virialization 
of the systems commences once K/\W\ ~ 0.5. (Bottom) Mi- 
gration paths followed by the aforementioned systems in a M 
vs Xmax plot. The paths converge to the line representing 
0-node equilibrium configurations, see Fig. I16l and Eq. 135H . 



will take a longer time for them to virialize completely. 

The examples with Vo = 0.93 (K/\W\ T=0 = 1.225) 
and 2p = 0.92 (K/\W\ T = = 1.252) illustrate how long 
it takes for a system to virialize as we consider lower 
values of i/jq. For instance, for ipQ = 0.92, the system 
will need much more time than shown in Fig. ll9l to begin 
to virialize, but we can anticipate that it will follow a 
similar path as ipo = 0.93. Actually, the migration paths 
of these systems were obtained by evolving them up to 
times of the order r > 10 4 . 

These examples show that all overcooled profiles 
(Jf/|W| < 0.5) rapidly settle down onto a 0-node equi- 
librium configuration, but overwarmed profiles require 
much longer times to stabilize, a fact that makes them 
difficult to evolve numerically. 

Despite of this, we conclude that if the overwarmed 
runs are followed during a sufficient long time, they will 
always find their way to an appropriate 0-node equilib- 
rium configuration. The reason for this is that there is 



an infinite number of 0-node equilibrium configurations, 
including those in which the scale parameter is too small, 
A« 1; that is, there is always an equilibrium configura- 
tion available for any evolved initial profile. 



VI. CONCLUSIONS 

In this paper, we have studied the evolution of a self- 
gravitating scalar field in the Newtonian regime using nu- 
merical methods. The latter were systematically tested 
for their accuracy, convergence and reliability. In all 
cases, the numerical code gave the expected correct re- 
sults of the 0-nodc equilibrium configurations and its 
perturbations. Also, the numerical code preserved the 
scaling invariance of the SN system all along the time of 
the evolutions. 

An important point was the boundary condition im- 
posed on the SN system. We found that the implemen- 
tation of a sponge by adding an imaginary potential is an 
appropriate boundary condition. It allowed us to main- 
tain under control the amount of scalar matter reflected 
by the numerical boundary. 

Our results imply that 0-node equilibrium configura- 
tions are intrinsically stable against radial perturbations, 
and that they play the important role of final states arbi- 
trary scalar systems settle down onto in the Newtonian 
regime. On the contrary, excited equilibrium configu- 
rations were found to be intrinsically unstable configu- 
rations. Even though they are initially virialized, they 
evolve towards a 0-node solution. 

An important point here is that, due to the scaling 
properties of the SN system, the study of the whole space 
of possible equilibrium configurations was reduced to the 
analysis of some properly sized configurations. Moreover, 
the same was done to study non-equilibrium configura- 
tions to give simplified and accurate simulations. 

Last but not least, we retook the analysis of the grav- 
itational cooling within the Newtonian regime of scalar 
fields. The main result is that, in the weak field limit, the 
gravitational cooling is a very efficient mechanism, which 
allows any initial configurations to decay into a 0-node 
Newtonian scalar soliton. So far, we have not found any 
evidence for systems that disperse away by ejecting all 
their mass. 

We expect that the results presented here would pro- 
vide useful information about structure formation in the 
universe, not only re garding models of scalar field dark 
matter as infH I5L l33l l34l |35|| . but also about other mod- 
els whose dynamics is governed by the SN system beyond 
spherical symmetry as inQ, case about which we expect 
to report in the near future. 
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